Least Squares Fitting-Perpendicular Offsets

Dimitris Sardelis-Theodoros Valahas

The aim of the present article is (i) to reproduce the results of the Mathworld presentation of the Least Squares Fitting-Perpendicular Offsets[1] by using a more standardized notation that simplifies matters considerably, (ii) to extend [1] by showing that only one of the two solutions derived there properly classifies as optimal-minimum, (iii) to complement [1] by deriving the appropriate coefficient of determination - the measure of the model's explanatory power- and the corresponding ANOVA for regression -the standard tool for performing hypotheses tests regarding the statistical significance of the relation between the variables involved, and, (iv) to compare the perpendicular offsets approach with the traditional approach mainly on efficiency grounds. In particular, a straight forward comparison of the determination coefficients of the two approaches reveals that the perpendicular offsets approach provides always a better fit. Finally, (v) a Mathematica coded demonstration follows which yields the essential parameters of the method for any given data.

Given any bivariate data one may define the following sums:

\(S\_xx ≡ ∑\+\(i\  = \ 1\)\%n\ \((x\_i - x\&_)\)\^2 = ∑\+\(i\  = \ 1\)\%n\ x\_i\%2 - n  x\&_\^2\)

\(S\_yy ≡ ∑\+\(i\  = \ 1\)\%n\ \((y\_i - y\&_)\)\^2 = ∑\+\(i\  = \ 1\)\%n\ y\_i\%2 - n  y\&_\^2\)

\(S\_xy ≡ ∑\+\(i\  = \ 1\)\%n\ \((x\_i - 
                x\&_)\) \((y\_i - 
                y\&_)\) = ∑\+\(i\  = \ 1\)\%n\ \(x\_i\) y\_i - 
        n \( x\&_\) y\&_\)

where \(TraditionalForm\`x\&_\) and \(TraditionalForm\`y\&_\) stand for the sample means for each variable. These sums are the numerators of the sample variances with respect to the variables x, y and of the sample covariance x - y.

The Error Sum of Squared Distances

The regression line being y (x) = a + b x, it is required to find a and b so that the sum of the squared distances of the data points from the regression line is a minimum. Following the standard statistical terminology, the latter sum should be called  "Sum of Squares for Error" and  it can be denoted as SSE. Since the distance of each data point ((Subscript[x, i]), (Subscript[y, i])) from the regression line is

\(d\_i = \(\(|\)\(y\_i - \((a + b\ x\_i)\)\)\(|\)\)\/\@\(1 + b\^2\)\)

SSE takes the form

\(SSE ≡ ∑\+\(i\  = \ 1\)\%n d\_i\%2 = \((1\/\(1 + b\^2\))\) \(∑\+\(i\  = \ 1\)\%n\((y\_i - \((a + b\ x\_i)\))\)\^2\)\)

which, upon using (1), (2) and (3), becomes

\(SSE = \((1\/\(1 + b\^2\))\) \((a\^2\ n + 2\ a\ b\ n\ x\&_ + 
            b\^2\ n\ x\&_\^2 - 2\ n\ a\ y\&_ - 2\ b\ n\ x\&_\ y\&_ + 
            n\ y\&_\^2 + b\^2\ S\_xx - 2\ b\ S\_xy + 
            S\_yy)\)\)

The First Partial Derivatives with respect to a and b

The first derivatives of SSE with respect a and b are

\(∂\_a 
        SSE = \(2\ a\ n + 2\ b\ x\&_\ n - 2\ y\&_\ n\)\/\(b\^2 + 1\)\)

and

\(∂\_b 
        SSE = \(2\ b\ S\_xx - 2\ S\_xy\)\/\(1 + b\^2\) - \(1\/\((1 + b\^2)\)\^2\) \((2\ b\ \((b\^2\ n\ x\&_\^2 - 2\ b\ n\ x\&_\ y\&_ + 
                  n\ y\&_\^2 + 2\ b\ n\ x\&_\ \((\(-b\)\ x\&_ + y\&_)\) - 
                  2\ n\ y\&_\ \((\(-b\)\ x\&_ + y\&_)\) + 
                  n\ \((\(-b\)\ x\&_ + y\&_)\)\^2 + b\^2\ S\_xx - 
                  2\ b\ S\_xy + S\_yy)\))\)\)

respectively. Setting  each of (7) and (8)  equal to zero and solving the system, one finds

\(\(b\_±\) = \(\((S\_yy - S\_xx)\) ± \@\(\((S\_yy - S\_xx)\)\^2 + 4\ S\_xy\%2\)\)\/\(2\ S\_xy\)\)

and

\(\(a\_±\) = 
      y\&_ - \(b\_±\) x\&_\)

Thus we have two relative extrema: [\(TraditionalForm\`\(a\_+\), \(b\_+\)\)] and [\(TraditionalForm\`\(a\_-\), \(b\_-\)\)] which correspond to two mutually perpendicular regression lines passing through the data centroid (\(TraditionalForm\`x\&_, y\&_\)). Calculating the corresponding errors, one finds

SSE [\(TraditionalForm\`\(a\_+\), \(b\_+\)\)] =(1/2) ((Subscript[S, xx]) + 
            (Subscript[S, yy]) - (Sqrt[4 (Subscript[S, xy]^2) + ((((Subscript[S, xx]) - (Subscript[S, yy])))^2)]))

SSE [\(TraditionalForm\`\(a\_-\), \(b\_-\)\)] =(1/2) ((Subscript[S, xx]) + 
            (Subscript[S, yy]) + (Sqrt[4 (Subscript[S, xy]^2) + ((((Subscript[S, xx]) - (Subscript[S, yy])))^2)]))

Hence, SSE[\(\(a\_+\), \(b\_+\)\)]≺SSE[\(\(a\_-\), \(b\_-\)\)]. Therefore, the solution [\(\(a\_+\), \(b\_+\)\)] must refer to a minimum.

The Second Derivatives with respect to a and b and the Test Discriminant

In view of showing that [\(\(a\_+\), \(b\_+\)\)] is indeed a minimum by applying the second derivative test, we  need to determine the second partial derivatives of SSE with respect to a and b as well as the corresponding test discriminant.

After a fair bit of algebra, the second partial derivatives are found to be

(∂_aa) SSE = ((2 n)/(1 + (b^2)))

\(∂\_ab 
        SSE = \(∂\_ba 
          SSE = \(2\ n\ x\&_\)\/\(1 + b\^2\) - \(2\ b\ \((2\ a\ n + 2\ b\ n\ x\&_ - 2\ n\ y\&_)\)\)\/\((1 + b\^2)\)\^2\)\)

\(∂\_bb 
        SSE = \(1\/\((1 + b\^2)\)\^3\) \((2\ \((a\^2\ \((\(-1\) + 
                      3\ b\^2)\)\ n + \((n - 3\ b\^2\ n)\)\ x\&_\^2 + 
                2\ b\ \((\(-3\) + b\^2)\)\ n\ x\&_\ \((a - y\&_)\) - 
                2\ a\ \((\(-1\) + 3\ b\^2)\)\ n\ y\&_ + \((\(-1\) + 
                      3\ b\^2)\)\ n\ y\&_\^2 + S\_xx - S\_yy + 
                b\ \((\(-3\)\ b\ S\_xx - 2\ \((\(-3\) + b\^2)\)\ S\_xy + 
                      3\ b\ S\_yy)\))\))\)\)

and, consequently, the test discriminant is

\(\(\(D ≡ \((∂\_aa SSE)\) \((∂\_bb 
                SSE)\) - \((∂\_ab SSE)\)\^2 = \(-\(1\/\((1 + b\^2)\)\^4\)\) \((4\ n\ \((a\^2\ \((1 + b\^2)\)\ n + 
                b\^2\ \((1 + b\^2)\)\ n\ x\&_\^2 + 
                2\ \((b + b\^3)\)\ n\ x\&_\ \((a - y\&_)\) - 
                2\ a\ \((1 + b\^2)\)\ n\ y\&_ + \((1 + b\^2)\)\ n\ y\&_\^2 - 
                S\_xx + S\_yy + 
                b\ \((3\ b\ S\_xx + 2\ \((\(-3\) + b\^2)\)\ S\_xy - 
                      3\ b\ S\_yy)\))\))\)\)\(
\)
    \)\)

The Second Derivative Test for [ \(a\_+\) ,\(b\_+\)]

For the + solutions (9) and (10), we have

(∂_a)SSE[ \(a\_+\) ,\(b\_+\)]=0

(∂_aa)SSE[ \(a\_+\) ,\(b\_+\)]=\(\(2  
                n\)\/\(1 + \((\(-S\_xx\) + S\_yy + \@\(4\ S\_xy\%2 + \((\(-S\_xx\) + S\_yy)\)\^2\))\)\^2\/\(4\ S\_xy\%2\)\)\)≻0

(∂_ab)SSE[ \(a\_+\) ,\(b\_+\)]= (∂_ba)SSE[ \(a\_+\) ,\(b\_+\)]=\(\((4\ n\ x\&_\ S\_xy\%2)\)/\((4\ S\_xy\%2 - \((S\_xx - 
                    S\_yy)\)\ \((\(-S\_xx\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\))\)\)

(∂_bb)SSE[ \(a\_+\) ,\(b\_+\)]=-\(\((8\ S\_xy\%2\ \((4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2)\)\ \((\(-S\_xy\%2\)\ \((\(-S\_xx\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\) + n\ x\&_\^2\ \((\(-2\)\ S\_xy\%2 + \((S\_xx - 
                              S\_yy)\)\ \((\(-S\_xx\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + 
                              S\_yy)\))\))\))\)/\((4\ S\_xy\%2 - \((S\_xx - S\_yy)\)\ \((\(-S\_xx\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\))\)\^3\)

Consequently, one finds that

D[\(a\_+\) ,\(b\_+\)]=\(-\(\((32\ n\ S\_xy\%6\ \((4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2)\)\ \((S\_xx - \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) - 
                  S\_yy)\))\)/\((\(-4\)\ S\_xy\%2 + \((S\_xx - S\_yy)\)\ \((\(-S\_xx\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\))\)\^4\)\)

Since (Subscript[S, xx]) - (Sqrt[4 (Subscript[S, xy]^2) + ((((Subscript[S, xx]) - (Subscript[S, yy])))^2)]) - 
          (Subscript[S, yy]) ≺ 0 for any (Subscript[S, xx])(≻0), (Subscript[S, yy])(≻0) and (Subscript[S, xy]), the discriminant in (21) is  positive definite. Therefore, the optimal solution [\(a\_+\) ,\(b\_+\)] indeed classifies as a relative minimum.

The Coefficient of Determination

In standard terminology, the total variation in the dependent variable of all data points from \(TraditionalForm\`y\&_\) is called "Sums of Squares Total" ( SST). In our case, SST denotes the total variation in both variables and it is defined in a natural way as the sum of the squared distances  of all data points from the centroid \(TraditionalForm\`\((x\&_, y\&_)\)\):

\(SST ≡ ∑\+\(i\  = \ 1\)\%n\ \((x\_i - x\&_)\)\^2 + ∑\+\(i\  = \ 1\)\%n\ \((y\_i - y\&_)\)\^2 = 
      S\_xx + S\_yy\)

Also following the standard terminology, one must define the "Sum of Squares for Regression" (SSR) as the variation of all data points accounted for by regression (explained variation), i.e., as SSR = SST - SSE. In our case,

\(SSR[\(a\_+\), \(b\_+\)]\)=(Subscript[S, xx])+(Subscript[S, yy])-SSE[\(\(a\_+\), \(b\_+\)\)]=(1/2) ((Subscript[S, xx]) + 
              (Subscript[S, yy]) + (Sqrt[4 (Subscript[S, xy]^2) + ((((Subscript[S, xx]) - (Subscript[S, yy])))^2)]))

Finally, the percentage of the total variation accounted for by regression, is tradinionally called  "coefficient of determination" and is denoted by (r^2), i.e., it is defined as (r^2) ≡ SSR/SST. In our case, the coefficient of determination is

\(r[\(a\_+\), \(b\_+\)]\^2 ≡ 
        SSR[\(a\_+\), \(b\_+\)]\/SST = \(\((S\_xx + S\_yy)\) + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\)\)\/\(2\ \((S\_xx + S\_yy)\)\)\)

ANOVA for Regression

To test whether or not the x and the y variables are significantly linearly related, it is common practice to use the ANOVA for regression scheme and compare the resulting F-ratio with a critical F-value depending on a pre-set (customarily large) degree of confidence. In our case, the ANOVA  table  for regression  has the form

Comparison with the Standard Regression Model

The plaussible reasons given in [1] for the almost universal adoption of the standard regression model are that (a) the model provides a fitting function for the independent variable X that is used to estimate y for a given x(most often what the experimenter wants), (b) it allows uncertainties of the data points along the x and the y axes to be incorporated simply, and (c) it provides a most  simple analytic form for the fitting parameters.

Despite the reasonable arguments  of simplicity, one cannot overlook the fact that the standard regression model discriminates the variables as dependent (y) and independent (x) in an arbitrary, observer-dependent,  subjective way that defies the logical and scientific content of the relation between the variables involved. Besides, there are many experimental occasions where both variables are simultaneously subject to error and consequently, they should be treated in a more unified, joint manner.

Finally, adopting the standard statistics practice of selecting among alternative models as best the one with the largest determination coefficient, one may compare the determination coefficient in (24)  with the one of the standard regression model., i.e., (r^2) = ((Subscript[S, xy])^2)/(Subscript[S, xx])  (Subscript[S, yy]). Thus, taking the ratio of the two, one finds

\(\(\(r\^2\/r[\(a\_+\), \(b\_+\)]\^2 = \((2\ S\_xy\%2\ \((S\_xx + 
                S\_yy)\)\ \((S\_xx - \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\))\)/\((S\_xx\ S\_yy\ \((S\_xx + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\)\^2)\)\)\(
\)
    \)\)

Since

\(S\_xx - \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + 
        S\_yy ≺ \((S\_xx + \@\(4\ S\_xy\%2 + \((S\_xx - S\_yy)\)\^2\) + S\_yy)\)\^2\)

and

(Subscript[S, xx]) + (Subscript[S, yy]) ≺ (Subscript[S, xx]) (Subscript[S, yy])

we have that

\(r\^2\/r[\(a\_+\), \(b\_+\)]\^2 ≺ 1\)

Thus, the perpendicular offsets model provides a better fit from the standard regression model.

A demonstration with Mathematica

What follows is intended to demonstrate the use of the method for computing the basic regression parameters using Mathematica.

data

As an example, let us take the following data:

data = {{0, 0}, {50, 7}, {70, 8.8}, {80, 9.9}, {20, 4}, {100, 
            11.9}, {180, 20}, {130, 14.8}, {160, 17}, {110, 
            12.6}};

definitions

Let the first elements in the lists be the x-values and the second elements be the y-values. Define

n = Length[data];
datax = data[[All, 1]];
datay = data[[All, 2]];
sumx = Apply[Plus, Part[datax]];
sumy = Apply[Plus, Part[datay]];
xbar = N[sumx/n];
ybar = N[sumy/n];
sumxy = Dot @@ Transpose[data];
sumxsq = Apply[Plus, Part[datax^2]];
sumysq = Apply[Plus, Part[datay^2]];
SXX = sumxsq - n*xbar^2;
SYY = sumysq - n*ybar^2;
SXY = sumxy - n*xbar*ybar;
b = 
        N[((SYY - SXX + (Sqrt[(SYY - SXX)^2 + 4*SXY^2]))/(2*SXY))];
a = ybar - b*xbar;
Rsquared = 
        N[(SXX + 
                SYY + (Sqrt[(SXX - SYY)^2 + 4*(SXY)^2]))/(2*(SXX + SYY))];
SSE = 
        N[(1/2) (SYY + 
                SXX - (Sqrt[(SYY - SXX)^2 + 4*SXY^2]))];
SST = SXX + SYY;
SSR = SST - SSE;
MSR = SSR/1;
MSE = SSE/(n - 2);

model parameters

Using the definitions and the data above, one then finds

TableForm[{{"n", n}, {"sumx", sumx}, {"sumy", 
          sumy}, {"sumxsq", sumxsq}, {"sumysq", 
          sumysq}, {"sumxy", sumxy}, {"a", a}, {"b", 
          b}, {"Rsquared", Rsquared}}]
\(TableForm[{{"ANOVA TABLE"}, {"Variation", "Sums", "df", "Means", "F-ratio"}, {"Model", SSR, 1, MSR, 
          MSR/MSE}, {"Error", SSE\ , Length[data] - 2, 
          MSE}, {"Total", SST, Length[data] - 1}}]\)

[1] Eric W. Weisstein. "Least Squares Fitting--Perpendicular Offsets." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/LeastSquaresFittingPerpendicularOffsets.html